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PROLOGUE 


The global oceans contain approximately 50 times as much CO 2 in dissolved forms as that in 
the atmosphere, while the land biosphere including the biota and soil carbon contains about 3 times 
as much carbon (as C0 2 ) as that in the atmosphere (Sundquist, 1985). Thus the spatial and temporal 
variability of C0 2 fluxes over the ocean are crucial for projecting the future level of atmospheric 

C0 2 . 

One of the most important parameters of the oceanic C0 2 system is the partial pressure of 
dissolved carbon dioxide, pC0 2 , in the surface ocean ( Takahashi et al., 1993; Wong and Chan, 

1991). The difference between pC0 2 in surface seawater and in the overlying atmosphere defines the 
source and sink areas of C0 2 over the global oceans. Since the high latitude waters are probably 
undersaturated with respect to C0 2 in the summer ( Keeling , 1968), these oceanic areas could play an 
important role in climate-C<9 2 feedback processes by removing large quantities of C0 2 from the 
atmosphere (Wong and Chan , 1991). Temperate and polar oceans of both hemispheres are the major 
sinks for atmospheric C0 2 , whereas the equatorial oceans are the major sources for C0 2 ( Takahashi 
et al., 1997). Thus, the evaluation of the air-sea exchange of C0 2 is crucial to determine local and 
global balances of carbon in the atmosphere-ocean system. 

The evaluation of the atmosphere-ocean C0 2 exchange is regulated by the gradient of pC0 2 
across the air-sea interface, the gas transfer velocity (or piston velocity), and the solubility of C0 2 in 
water. There are a few methods available for evaluating the transfer velocity of C0 2 air-sea ex- 
change, which can be obtained by field measurements ( Broecker and Peng, 1974) or in the labora- 
tory (Liss, 1988). Field methods were applied using a variety of data sets obtained from numerous 
experiments, e.g., Barbados Oceanographic and Meteorological Experiment (BOMEX), Geochemi- 
cal Ocean Sections Study (GEOSECS), and Transient Tracers in the Ocean (TTO) programs 
(Broecker and Peng, 1971; Peng et al., 1974; Peng et a\., 1979; Smethie et al., 1985, Batrakov et cil., 
1981). Laboratory methods using wind/wave tunnel experiments can be used to relate the wind 
speeds with the measured transfer velocities (Liss and Merlivat, 1986). This was the approach 
employed in the Programme Fran$ais Ocean-Climat en Atlantique Equatorial (FOCAL) described by 
Andrie et al. [1986], Oudot andAndrie [1986], and Oudot et al. [1987]. In addition to the 
GEOSECS and TTO programs, there are also some studies on the variability of the oceanic C0 2 
system in the subarctic Atlantic Ocean ( Takahashi et al., 1983', Takahashi et al., 1985), the tropical 
Atlantic Ocean (Oudot andAndrie 1989), the western Pacific Ocean (Fushimi, 1987', Inoue et al., 
1987), the Southern Ocean (Inoue and Sugimura, 1988), and the subarctic North Pacific (Gordon et 
al, 1971; 1973; Takahashi, 1989; Murphy et al, 1998). 

The seasonal and interannual variations of C0 2 in the surface oceans are not only affected by 
the air-sea exchange physical processes but also by the photosynthetic uptake of C0 2 by phytoplank- 
ton. For instance, spring phytoplankton blooms in the surface waters of the North Atlantic Ocean 
can cause a precipitous reduction of surface water pC0 2 , C0 2 and nutrients in a span of a couple of 
weeks. The mechanisms that drive this large biogeochemical variability were modeled by previous 
investigators (for example; Taylor et al, 1991). In contrast, seasonal changes in C0 2 and nutrients 
are more gradual in the North Pacific and macro-nutrients are only partially consumed in the surface 
waters of the subarctic North Pacific Ocean (Takahashi et al., 1993). 
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In this TM, we describe the model functionality and analyze its application to the seasonal and 
interannual variations of phytoplankton, nutrients, pC0 2 , and C0 2 concentrations in the eastern 
subarctic Pacific at Ocean Weather Station P (OWS P, 50° N 145° W). We use a verified one-dimen- 
sional ecosystem model ( McClain et al., 1996), coupled with newly incorporated carbon flux and 
carbon chemistry components, to simulate 22 years (1958-1980) of pC0 2 and C0 2 variability at 
Ocean Weather Station P (OWS P). This relatively long period of simulation verifies and extends the 
findings of previous studies (Wong and Chan, 1991 ; Archer et al., 1993; Antoine and Morel, 1995a; 
Antoine and Morel, 1995b) using an explicit approach for the biological component and realistic 
coupling with the carbon flux dynamics. The slow currents and the horizontally homogeneous ocean 
in the subarctic Pacific make OWS P one of the best available candidates for modeling the chemistry 
of the upper ocean in one dimension. The chlorophyll and ocean currents composite for 1998 shown 
in Figure 1 illustrates this premise. The chlorophyll concentration map was derived from SeaWiFS 
data and the currents are from an OGCM simulation (R. Murtugudde, personal communication). 
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Figure 1. Chlorophyll and ocean currents composite for 1998 based on SeaWiFS data and 
OGCM simulations, respectively. The black dot denotes the location of OWS P. 
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1.0 MODEL DESCRIPTION 


The model consists of a previously verified ( McClain et al., 1996) four-component ecosystem 
model, coupled with newly added carbon dioxide and oxygen components. The ecosystem model is 
physically forced by sea surface temperature (SST) and mixed layer depth values originating from the 
Garwood mixed layer model {Garwood, 1977). A description of the mixed layer model validation and 
application to the forcing of the ecosystem model is given in McClain et al. [1996] and is not repro- 
duced here. The ecosystem model has some modifications which are explicitly defined in McClain et 
al. [1998] and McClain et al. [1999]. In addition to these modifications, two new components, iron 
(Fe) and dissolved organic carbon (DOC), were added to the model. Previous studies {Maldonado et 
al., 1999; Boyd et al, 1996; Martin et al, 1989 ) showed that the phytoplankton growth in the subarctic 
Pacific is generally iron limited. We added the iron component to obtain more accurate simulations of 
nitrate, phytoplankton and community production. The DOC component was added to more accurately 
reproduce the seasonality of the net community production via the respiration of DOC by bacteria 
during the fall-winter period. In addition, the light penetration model uses a different spectral set of 
water absorption (a w ) coefficients and phytoplankton specific absorption (a* p ) coefficients. Figure 2 
shows the spectral variability of these coefficients. The clear sky irradiance is modified to account for 
the observed cloud cover by applying a power law correction ( McClain et al., 1996) tuned to yield the 
observed climatological monthly mean surface irradiances {Dobson and Smith, 1988). The climato- 
logical monthly irradiance, obtained by averaging the daily mean irradiance for each month and for the 
entire period of simulation (1960- 1 980), is shown in Figure 3. The seasonal variability shown in Figure 
3 agrees with the values given by Dobson and Smith [1988]. The coefficient for the cloud attenuation 
formulation of McClain et al. [1996], given by 

E = E^M^Clcn (1) 

was changed from 0.53 to 0.45 to better match the climatological cloudy sky radiation {E) of Dobson 
and Smith [1988]. This new coefficent also provided an improved match of the simulated SST with the 
observed values. The cloud cover {Cld, in oktas) was obtained from OWS P observations. 


The following system of coupled differential equations simulates the dynamics of phytoplankton 
nitrogen {P), zooplankton nitrogen (Z), ammonium {NH 4 ), nitrate (NO s ), iron {Fe), dissolved organic 
carbon {DOC), total carbon dioxide {C0 2 ), and oxygen {0 2 ) stocks within the upper ocean: 
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Figure 2. Spectral light absorption coefficient for sea water (a w (A)) and specific chlorophyll ab- 
sorption coefficient (a* (X)) adopted in the spectral downwelling irradiance model formulation. 



Figure 3. Comparison between the climatological daily averaged cloudy sky surface solar 
irradiance from Dobson and Smith [1988] and model simulation. 


2 




( 6 ) 


dFe ... dFe d 
dt ' dz dz 


dFe 

dz 


=S(z)FF e -N p % 


dDOC , dDOC d 

~-dT +w '-dT'dz 


K 


dDOC. 
dz . 


= ( a'pmP+a' z (g+hZ)Z)jj - k rc DOC 


dTC0 2 

~ aT 


dTCO 


■+w. 


2 d 


dz dz 


K- 


dTCO* 


dz 


=5(z)FC0 2 -N p ±+k rc D0C 


do 2 

~df 


do 


' 2 d 


■ +w.^— - -z— 


dz dz L ’ dz J 


K, 


dOi 1 


=5(z) F0 2 + 


N PT k rc D0C 


O' 


(7) 

( 8 ) 

(9) 


A flow chart of the coupled model, showing the interaction between its major components, is 
shown in Figure 4. Table 1 summarizes the parameters used in Equations (2) through (9) and derived 
quantities. The ratios C/N, 0 2 A C are the carbon-to-nitrogen and oxygen-to-carbon Redfield ratios ( 106: 1 6, 
138:106), respectively. The iron-to-nitrogen ratio is assumed to be 19.8 |imol/mol. The Kroenecker 
delta (5[z=0]=l; 5[z>0]=0) is used to denote that the carbon dioxide and oxygen fluxes ( FC0 2 and 
FO , respectively) are applied at the air-sea interface only. Details of the Ekman upwelling (w) and 
vertical diffusion (X) formulations, the numerical method to solve the coupled differential equations, 
and the various parameters used by the ecosystem model components are given in McClain et al. 
[1996]. The Ekman upwelling profile function was modified to account for depth attenuation below 
the Ekman depth ( D e ). The formulation for D e ( Pond and Pickard, 1978) as a function of wind speed W 
and latitude/, the shape function R(z), and the vertical velocity are given by: 



Figure 4. Flowchart showing the principal components of the ecosystem/carbon-flux 

one-dimensional model. 


3 




(10a) 


'.-(f) 


1/3 


D. 


D - 


/ 

4JW 


' [sin(H)]° 5 

-7TZ 

R(z)= 1 -cos [ ) e ° e > for z < D t 


f z-D' V 5 1 

k(z-D') ' 

(z-D COS 

v max e / -1 

L2(z m „- D')\ 


, for z> D e 


(10b) 

( 11 ) 

(12) 
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Equation 10b (Pond and Pickard, 1978) produces nearly identical values of D e as the standard 
Ekman formula given in Equation 10a. The standard formula introduced numerical instabilities, so 
Equation 10b was used. The shape function R(z) is constructed based on McClain et al. [1996] for the 
portion above D t (equation 11). For z > D e , R(z) is formulated such that continuity at the inflection 
points z=D e and z-z.^ is maintained, and that R(z) decreases to zero monotonically towards z=z max = 350 
meters (equation 12). 


The net community production, N , is defined as the summation of all sources and sinks of nitrate 


and ammonium in equations (3) and (4): 
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The ratio N/C is 16:106 (Redfield ratio). The effective growth is a function of the maximum 
growth (G o ), light limitation (L //m ), and nutrient limitation (N lim ): 


r a r & gP ^ 

G-fi G 0 e SP 

£,im='-“p( 


r e . i 

in 

\ 

V max 


Fe 


Fe \im ( K F +Fe) 

N y =NHy +NO t y 
lim mm 3hm 


NO^.= 


Nt > 


(1 -NH 4 ) 


3lim " ^no3 + ^3> 


NH = 


NH. 


4lim ~ (K Nm +NHf) 


(15) 

(16) 

(17) 

(18) 

(19) 

( 20 ) 
( 21 ) 


4 


2.0 MODEL FORCING AND BOUNDARY CONDITIONS 


The Neumann boundary condition, dXJd z=0, is applied at both the surface and lower (350 m) 
model domain boundaries for P and Z. Initial profiles of temperature and N0 3 were obtained from 
winter and annual climatologies ( Conkright , et al, 1994), respectively. Depth-independent initial con- 
centrations of P and Z are 0.2 and 0.1 mmol N nr 3 , respectively. Observed fall profiles of NH 4 were 
used as initial conditions (Frost, 1993). At the lower boundary, fixed values equal to 0.0 and 41 .5 mmol 
N nr 3 were applied to NH 4 and NO r respectively. The Neumann condition was applied at the surface. 

For TCO and 0 2 , climatological values at the surface and at the bottom were used to construct 
the initial linear profiles (surface TCO =2050 mmol nr 3 and bottom C0 2 =2100 mmol nr 3 ; surface 
O =320 mmol nr 3 and bottom O 2 =100 mmol nr 3 ). A fixed value equal to the initial condition was 
applied at the lower boundary. The following formulations for the C0 2 and 0 2 gas exchanges were 
applied in the form of flux boundary conditions ( FC0 2 and F0 2 in mmol nr 2 d 1 ) at the sea-air interface: 


FCO=K o a&pC0 2 ( 22 ) 

F0=K o [0 2 (sat}-0 2 ] (23) 

a = exp[-60.2409+ 9345.17 + 23.35851og + 5(0.0235 17-0.0236567’+0.0047036T) 2 ] (24) 


where, k is the gas exchange coefficient (piston velocity, in m/d) which is a function of water tempera- 
ture and wind speed (Wanninkhof, 1992; Liss and Merlivat, 1986), a is the C0 2 solubility in seawater 
(in mmol m' 3 (latm -1 ) which is a function of temperature and salinity (Weiss, 1974), ApC0 2 (in |iatm) is 
the difference between air and sea pC0 2 , and 0 2 (sat) is the oxygen saturation concentration (in mmol 
nr 3 ) in seawater which is a function of temperature and atmospheric pressure (Weiss, 1970). 

The eolian iron flux, FFe, is given by: 

C . Sc F Pr Cp Sp 

J7T7p — dus< Fe O Is. (25) 

55.84 

where, C is the seasonal atmospheric dust concentration in ng/(kg of air) from Tegen and Fung 
[1994], Sc F = 1000 is the scavenging ratio (Duce, 1995), Pr is the daily National Centers for Environ- 
mental Prediction (NCEP) precipitation in kg/m 2 /d, C Fe =0.035 is the iron mass fraction in the dust, and 
S =0.1 is the soluble iron fraction. The numerical factors are used to convert dust concentration to an 
iron flux (FFe) in pmol/m 2 /d. The seasonal variability of precipitation and the iron flux for the period 
of 1970-1980 is shown in Figure 5. 
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Table 1. Summary of model variables and input parameter definitions and values. Initial surface and 

bottom concentrations are provided when appropriate 



Fe 


DOC 


Chl-a/N 


Value 


0.2 


0.1 


0.1 


Climatology 


50/400 


10/0 


2050/2100 


320/100 



0.05 


0.5899 


0.0633 


0.02 


0.0633 


50 (<60m) 



0.3 


0.019 


Definition 


Initial phytoplankton concentration (mmol N/m 3 


Initial zooplankton concentration (mmol N/m 3 


Initial ammonium concentration (mmol N/m 3 ) 


Initial nitrate concentration (mmol N/m 3 ) 


Initial iron concentration (pM' 


Initial dissolved organic carbon (mmol C/m 3 


Initial total CO, Concentration (mmol C0 2 /m 3 ) 


Initial oxygen concentration (mmol OJ m 3 ) 


Regenerated production fraction 


New production fraction 


Half saturation for NO, uptake (mmol N/m 3 ) 


Half saturation for NH 4 uptake (mmol N/m 3 ) 


Half saturation for Fe uptake (pM) 


Maximum zooplankton grazing rate (a 


Net community production (mmol C/ m 3 /d) 


Phytoplankton death rate (d *) 


Phytoplankton growth rate at 0 °C (d ‘) 


Temperature sensitivity of algal growth (°C ') 


Respiration coefficient 


Temperature sensitivity of algal resp. (°C 1 ) 


Maximum photoacclimation param. (uEin/m 2 /s) 


Maximum photoacclimation param. (pEin/m 2 /s' 


Ammonium inhibition of NO s uptake 


Maximum phytoplankton sinking Rate (m/d) 


Chlorophyll to nitrogen ratio 


Zooplankton death rate (d ‘ 


uadratic coefficient of zooplankton mortality (mmol' 1 m 3 d' 1 


Maximum zooplankton grazing rate (d -1 


constant (m 3 /mmol N 


nilB5IHSIiSW5Blllll.MlIJJ.rai 



ItSEMUnuitli 


Unassimilated zooplankton ingestion Ratio 


Respiration rate for zooplankton at 0 °C (d -1 


Temperature sensitivity of Z respiration (d r 


Fraction of dead P converted to NH 4 


Fraction of dead Z converted to NH 4 


Fraction of dead P converted to DOC 


Fraction of dead Z converted to DOC 


Respiration rate of DOC into CO 2 (d 


Maximum rate of nitrification (mmol/d 


Minimum inhibition dosage for nitrification (W/m 2 


Half saturation dosage for nitrification photoinhibition 


Fecal pellet remineralization fraction 


Minimum bottom eddy diffusion (m 2 /d) 
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NCEP (50.5°N,144.4°W) 



Figure 5.Time series of precipitation rate and eolian iron flux for 1970-1980. 


The 30-year pC0 2 atmospheric time series used to force the model is based on a least-squares fit to 
the 15 years (1978-1993) of pC0 2 observations in Cold Bay, Alaska, the closest long-term monitoring 
site at roughly the same latitude of OWS P. The pC0 2 time series is given by: 

pC0 2 (air)=A o +A x t+A 2 sm +A 3 ) +A 4 sin +A 5 ]+A 6 sin py +A ? ] (26) 

where t is the time in months, A 0 = 280.8 patm is the intercept, A = 0.134 patm mo' 1 is the slope, the 
amplitudes ( A 2 , A 4 , A 6 ) are 6.61, 3.01, and 0.87 patm, respectively. The phases (A 3 , A 5 , A 7 ) are 0.77, - 
15.2, and 0.15 radians, respectively. Figure 6 shows the time series of atmospheric pC0 2 at Cold Bay, 
at Mauna Loa, and the synthesized signal (dashed line) using the analytical formula (26). Note that the 
long-term trend at Cold Bay and Mauna Loa are essentially identical, while there is a significant differ- 
ence in the seasonal amplitudes. The fact that the trends are similar and the amplitudes are different can 
be attributed to atmospheric mixing processes. Specifically, the atmospheric seasonal mixing is im- 
parted preferentially along latitude lines by virtue of stronger zonal flows in the atmosphere, allowing 
meridional gradients of pC0 2 to be much stronger than the zonal gradients. Conversely, the long-term 
trend in pC0 2 distribution is latitudinally more uniform because the meridional mixing time scale is 
short when compared to the pC0 2 trend due to antropogenic sources. 

Sensitivity tests were conducted off-line using the air-sea flux formulation implemented in the 
model. Four formulations were tested; Liss and Merlivat [1986], Tans etal. [1990], Wann inkhof [1992], 
and McGillis and Wanninkhof [1999]. Monthly averaged (1973-1978) winds, SST, salinity, and ApC0 2 
from Wong and Chan [1991] were used in these calculations. The Liss and Merlivat [1986] method 
provides the smallest sea-air flux (-2.19 mmol C0 2 nr 2 d 1 ), whereas the methods of Tans et al. [1990] 
and Wanninkhof [1992] provide very similar results (Figure 7). The largest sea-air flux (-31.8 mmol 
C0 2 m' 2 d 1 ) was provided by the McGillis and Wanninkof[ 1999] method since it is a cubic function of 

the wind speed (W): 
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Months after January 1, 1978 


Figure 6. Time series of Mauna Loa and Cold Bay atmospheric pC0 2 . 
The dashed line represents the time series of pC0 2 used to force the model. 

Comparison of Gas Exchange Coefficient Formulations 
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ApC0 2 from Wong and Chan, 1991 (patm) - Solid 
Climatological Winds (1958-1980) - Dotted 



1) Sea-Air Flux (-2.19 mmolCOj/mVd) from Liss and Merlivat [1986] 

2) Sea-Air Flux (-4.54 mmoICCyn^/d) from Tans et al. [1990] 

3) Sea-Air Flux (-4.76 mmoIC(ym 2 /d) from Wan ninkhof [1992] 

4) Sea- Air Flux (-31.77 mmolCOj/m^d) from Wannlnkhof and McGIllis [1999] 


Figure 7. Comparison between sea-air C0 2 flux calculations using four different 
parameterizations for the gas exchange coefficient. 




( 27 ) 


( _ v-I/Z. 
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where Sc is the Schmidt number given as a function of temperature T as: 

Sc = 2073. 1- 125. 6274-3. 627T 2 -0.043219r 3 (28) 

Equation (27) is used for the stand-alone flux calculations using climatological winds and various 
formulations for K. For the actual model simulations we use the formula suggested by Wanninkhof and 
McGillis [1999] for short-term (< 1 day) winds given by: 

. . - 1/2 

K=0.0283lV 3 (^) (29) 

We adopted the McGillis and Wanninkhof, [1999] method for this study. We conducted several 
sensitivity runs using different methods to derive the gas exchange coefficient. Figure 8 shows the 
ApCO and air-sea flux seasonal plots simulated by the model for three different gas exchange coeffi- 
cient methods. The surface C0 2 fluxes obtained with the Liss and Merlivat [1986], Wanninkhof [ 1992], 
and Wanninkhof and McGillis [1999] are -6.8, -8.8, and -10.5 mmol/m7d, respectively. The negative 
sign indicates that the flux is from the atmosphere to the ocean in all three cases. We discuss more 
details of these results in section 4. 




1) Sea-Air Flux (-6.80 mmoICCyitf/d) from Llss and Merlivat [1986] 

2) Sea- Air Flux (-8.75 mmo\COJm 2 fd) from Wanninkhof [1992] 

3) Sea-Air Flux (-10.51 mmolCO^m^d) from Wanninkhof and McGillis [1999] 

Figure 8. Seasonal variability of sea-air &pC0 2 and surface flux predicted by the model using 

three different gas exchange formulations. 
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3.0 OCEANIC pC0 2 FORMULATION 


To calculate the pC0 2 concentration in seawater we must first understand the thermodynamics of 
the C0 2 system. The total C0 2 concentration in seawater, TCO r can be written as: 


tco 2 



+ 



(30) 


where all quantities in square brackets are stoichiometric concentrations; [CC> 2 ] is the dissolved carbon 
dioxide, [. HCO ’ *] is the bicarbonate ion, and [CO/] is the carbonate ion. Another quantity that influ- 
ences the calculation of pC0 2 is the total alkalinity, TA, which can be written as the sum of its major 
terms: 


TA =CA+ 


B{OH) a 


+ 


OH' 



(31) 


where CA=[HCO s ]+2[CO /'] is the carbonate alkalinity, [B(OH) 4 ] is the borate ion concentration, and 
[OH ] and [ff + ] are the products of H 2 0 dissociation. The partial pressure of dissolved C0 2 is defined 
by the relationship: 


pC0 2 = 


[C0 2 ] 

a 


(32) 


where [CO/ represents the carbon dioxide in solution and a (mmol nr 3 jiatm 1 ) is the solubility of 
carbon dioxide in seawater (Weiss, 1974). We can calculate [C0 2 ] from TC0 2 and TA (the values for the 
hydrogen ion activity, H + , and the carbonate alkalinity, CA, are also computed). We follow the recur- 
sive formulation of Antoine and Morel [1995] to estimate the dissolved [C0 2 ] concentration 


[C0 2 ]=TC0 2 -A+ 


fAK r -TC0 2 K r -4A+Zl 

2(K r -4) 


(33) 


where A takes the value of TA initially. K r is equal to K } /K 2 , where K t and K 2 (mmol nr 3 ) are the 
dissociation constants of carbonic acid, which are a function of temperature and salinity ( Goyet and 
Poisson , 1989). The other quantities required to calculate [C0 2 ] are 

Z = 4 (TC0 2 Kf+AK r (2TC0 2 - A)(4-K r ) (34) 


+ _([C0 2 ]K { ) laCOJKf+SAlCOJKfo 
~ 24 24 


CA=TA-BR 


h 

K b+ H+ 


H + 


(36) 


K x =I0' pK \K 2 =10' pK2 

pK x = m T 21 - + 3.356 - 0.00171 S log(T) + 0.000091 S 2 
pKi= M50J*? + 4 . 6 04 . 0.00385 S log (T) + 0.000182 S 2 


(37) 

(38) 

(39) 
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(40) 


K. 


—10 


pK w = M + 2.241 - 0.09415 S 0 5 


(41) 


, -8966.9 - 28 90.53S 05 - 77.942S+1.728S 15 - 0.099635 2 . 

— j ^ 

1 48.0248+ 137.1 9425 0 5 + 1 .62 1 425+ 

(-24.4344-25.085S a5 -0.2474S)log(7)+0.0531057iS' 0 ' 5 (42) 

where BR = 0.00042S/35.0 mmol nr 3 is the borate concentration ( Whitfield and Turner, 1981), and K B 
and K w are the dissociation constants for boric acid and sea water, respectively (Weiss, 1974). 

The above equations are repeated, with the variable A taking the value of CA at the previous 
iteration, until | CA - A\ < 10 8 . The value of pC0 2 is then estimated according to: 


pC0 2 = 


TC0 2 


a 



Mi 

H +2 


(43) 


Table 2 shows a comparison between pC0 2 values measured at OWS P (Wong and Chan, 1991) 
and those calculated using equations (33) through (44). The data were acquired during the winter of 
1975 (17 January and 9 February). A very close agreement (less than 1%) is achieved between the 
observed and calculated pC0 2 values. 


Table 2. Oceanic properties of surface sea water from OWS P in early 1975 (Wong and Chan, 1991). 
The measured (a) pC02 values are compared with those calculated (b) using equations (30) through 
(43). The difference between measured and calculated values (e) is 2.1 patm 


Date 

SST 

SSS 

TA 

TC02 

rrrnMi 

P co 2 b 

£ | 

17 Jan 

5.8 

32.656 

2216 

2039 


310.5 

-2.1 


5.8 

32.678 

2203 

2030 

317.4 

315.3 

-2.1 


Data collected during 1973-1978 (Wong and Chan, 1991) were used to validate and formulate 
model parameters. The salinity data is much sparser so that an analytical formulation must be obtained 
to provide an uninterrupted sequence of values to force the model. A regression of salinity versus 
temperature yields the linear relationship S = 32.8124 - 0.017197’. The total alkalinity can be expressed 
as a linear function of salinity. For example, the linear formulation for TA of Clayton et al. [1995], with 
an intercept adapted to reproduce the mean TA values at OWS P, can be used to produce an hourly time 
series of TA values to force the model (TA = 1 17.8 + 64.2325). However, since the salinity at OWS P 
does not change significantly (32.59 to 32.72 psu), the model results are not significantly different than 
assuming a constant value for TA. Also, there is little information available on the seasonal variability 
of TA at OWS P for model verification. Therefore, we used a constant value of 2215 |ieq/kg, which is 
also the value used by Antoine and Morel [1995b], and calculated at OWS P in January 1975 by Wong 
and Chan [1991]. Table 3 summarizes the parameters used in the geochemical portions of the model 
(iron flux and carbonate chemistry). 
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Table 3. Summary of parameters used in the iron flux formulation and carbonate chemistry 



Value 

Definition 


Climatological 

Atmospheric dust concentration Qig/kg of air) 

Sc Fe 

1000 

Iron scavenging ratio 

C F ,. 

0.035 

Iron mass fraction 

s Fe 

0.1 

Soluble iron fraction 

Ko 

Eqs. 27, 29 

Gas exchange coefficient 

Sc 


Schmidt number 

TA 

2215 

Total Alkalinity (|ieq/kg) 

B339H 

Eq. 33 

Carbon dioxide in solution (mmol/kg) 

H 

Eqs. 32, 43 

Partial pressure of C0 2 in sea water (jiatm) 

a 

Eq. 24 

CO 2 solubility in sea water (mmo!/rn 3 /jiatm) 

CA 

nmm mmm 


cmm 


Dissociation constants for carbonic acid 

UL« 

Eqs. 40, 42 

Dissociation constants for borate and water 

BR 

0.00042S/35 

Borate concentration (mmol/m 3 ) 

H + 


Hydrogen ion concentration 


4.0 MODEL SKILL ASSESSMENT 

4.1 Model Sensitivity to Forcing and Parameterization 

Table 4 shows a summary of the model response to different forcing and parameterization. We 
provide a comparison between the model results and available data sets from OWS P in the next sub- 
section. An additional parameter test was performed to evaluate the sensitivity of the carbon dioxide 
variables to changes in salinity. One of the objectives of this test is to evaluate the required sensitivity 
and accuracy of salinity remote sensors so that measurements can be useful in the determination of the 
surface carbon flux in the ocean. 


Table 4. Summary of model sensitivity to different parameterizations. The criteria for improvements 
include agreement with data such as phytoplankton (P) concentration, primary productivity 

(PP), and pC0 2 


Parameter 

Range 

Result 

Rm 

3.0-5.0 

4,0 for best P, PP 

a 

£? 

Linear to quadratic 

g=0.35Z 2 for best P, PP 

Fe, DOC 

Added components 

Improved PP seasonality 

a „(k) and a* p (X) 

For OWS P and Warm Pool 

Values for OWS P were used 

EJX) 

Up to 14% adjustment 

Improvement in SST 

Gas exchange K^W 2 ) 

1, 2, and 4xWanninkhof 

4x best pC0 2 variability 

K(w 3 ) 

K o =(0.0283W 3 )(Sc/660)‘°- 5 

Best agreement with 
observed pC0 2 

pC0 2 calculation from 
SST, S, TA and TCO } 

Antoine and Morel (1995), 
P.Murphy (personal 
communication), DOE 

All within less than 5% 

TC0 2 (350m) 

Constant to seasonal 

TC0 2 = 2100 best profile 

C0 2 flux (FCOj) forcing 

Top layer, distributed 
within mixed laver (ML) _ „ 

ML approach required for 
stabilitv (larse Ko onlv) 

Iron flux (FFe) forcing 

S F .=0.05,0. 1,0.2, 0.3 

Best P, PP with S Fe =0.05 

Carbon Production 

Cloern et al. (1995), Redfield Ratio 

Redfield Ratio more 
internally consistent 
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We use the relationship of Millero et al. [1998], which depends on temperature ( T) and salinity 
(5), to calculate the total alkalinity (7A) at each model time step 

NTA = 2300-7.00(7 -20)-0. 158(7-20) 2 
TA=NTA-^ 

where NTA is the salinity normalized alkalinity. The salinity is calculated from the regression on tem- 
perature based on the salinity and temperature data from OWS P 

5=32.8124-0.017197 (46) 

We conducted two runs, one using the salinity derived from equation (46), and another using the 
salinity in equation (46) with a positive bias of 0.5 psu. We adopted the gas exchange coefficient of 
Wanninkhof and McGillis [1999] in both runs. The results are summarized in Table 5. The parameters 
that are most sensitive to the salinity variability are the ApC0 2 , and, consequently, the air-sea C0 2 flux. 
These are also the two variables with standard deviations larger than the mean. The dissolved C0 2 does 
not change much (< 1%). The 24.3% increase in C0 2 uptake from the atmosphere, due to the 0.5 psu 
increase in salinity, is compensated by an approximately equivalent decrease in the C0 2 supply via the 
bottom flux. 

Note that the climatological ApC0 2 mean is reduced with the increased salinity, but the surface 
air-sea flux into the ocean actually increases despite the very small change in the gas exchange coeffi- 
cient. This result seems counter-intuitive at first glance but becomes clear after inspection of the sea- 
sonal variations of ApC0 2 and K. Figure 9 shows the seasonal variability of the surface C0 2 flux, 
ApCO and K. The solid and dashed lines represent results from runs 1 and 2, respectively. The K 
values 2 for the two runs are virtually identical. Figure 9 clearly shows that the ApC0 2 is largely positive 
from May through September and negative during the rest of the year. However, K is largest during the 
winter, spring and fall vigorous wind forcing, but significantly reduced during May through August 
when the wind speeds subside at OWS P. This phase opposition of ApC0 2 and K explains why the 
surface flux increases despite a climatological decrease in the ApC0 2 mean value. Note also that the 
ApC0 2 from Run 2 is more negative during winter, spring and fall than the ApC0 2 from Run 1. 



Table 5. Summary of sensitivity run based on a 0.5 psu salinity increase from Run 1 to Run 2. 

The averages are based on years 1960-1980 



Run 1 

Run 2 

Difference 

Percent 

S(psu) 

32.67 

33.17 

0.5 

1.5 

rCO,(|imol/kg) 

2007.1 ±22.7 

2033.1 ±22.6 

26.0 

1.3 

pC0 7 (\i atm) 

315.6 ±23.5 

314.2 ±24.4 

-1.4 

0.5 

FCO,(gC/m 2 /yr) 

-42.4 ± 47.9 

-52.7 ±55.6 ! 

-10.3 

24.3 

Bottom C0 2 Flux (gC/rrr/yT) 

-25.1 ±8.2 

-15.8 ±5.5 

9.31 ~ 1 

37.0 

TACueq/ks) 

2201.2 ±11.4 

2234.9 ±11.5 

33.69 

1.5 

CChi \ imnl/kcr^ 

14.8 ±0.9 

14.7 ±0.9 



Jil 

Ci( mmnl/kp/afrnl 

47.13 ±3.98 

47.0 ±4.00 

-0.13 

JL2 

AnCO,( uatrrO 

3.58 ±26.77 

2.15 ±27.75 

-1.43 

39.9 

tf(mmol/m 2 /d/^atm) 

1.376 ±2.752 

1.372 ±2.744 

-0.004 

0.3 
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Figure 9. Seasonal variability surface flux, sea-air ApC0 2 , and gas exchange coefficient 
predicted by the model for two runs with a salinity difference of 0.5 psu. 






4.2 Model-Data Comparison 


Table 6 provides a comparison between model and observed (Wong and Chan, 1991) biogeochemi- 
cal parameters. Both model and observed values are averaged over years 1973-1978. Primary produc- 
tion, chlorophyll and nitrate are depth-averaged (0-80 m). This particular model run uses the gas ex- 
change coefficient (K o ) from Wanninkhof [1992], The model primary production is within 20% of the 
observed value. Chlorophyll, nitrate, and SST are all within less than 5% of the observed values. The 
salinity is obtained via a regression with SST; thus it is not independent of the observations and should 
not be evaluated as a model parameter. A comparison between model (solid line) and observed (dotted 
line) seasonal temperature profiles, averaged over the overlapping model-data years (1958-1966), is 
shown in Figure 10. The agreement is quite good. The model profiles represent the seasonal variability 
of the temperature stratification and SST at OWS P reasonably well. A very small warming trend 
(0.001 °C/month not shown here) was identified below 175 meters (below the euphotic zone). We 
removed the trend after the fact for display purposes only. These small discrepancies do not signifi- 
cantly affect the model results. 

Winter T(°C) Spring T (°C) 

23456789 10 11 12 23456789 10 11 12 



Figure 10. Comparison between model and observed seasonal temperature profiles 

averaged over the period of 1958-1966. 
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The model TC0 2 is about 2% lower on average than the observed TCO r The mean atmospheric 
pC0 2 is within 3% of the measured mean. The largest uncertainty lies in the calculation of the air-sea 
C0 2 flux ( FC0 2 ). The surface flux is the most elusive quantity in the oceanic C0 2 measurements and 
modeling. It depends on accurate measurements of ApC0 2 and a robust formulation for the gas ex- 
change coefficient (K g ), since direct measurements of the air-sea flux are not yet available. As we have 
shown in section 2, there are multiple gas exchange coefficient formulations in the literature which 
yield a wide range of values for a given wind speed. For example, as shown in Table 7, FC0 2 can 
change significantly depending on which formulation for K o is used. The value calculated from the 
observed ApC0 2 and K g from Liss and Merlivat [1986] is -8.3 gC/m 2 /yr; if K g from Wanninkhof [ 1 992] 
is used the flux increases to -19.1 gC/m 2 /yr. Alternatively, the flux calculated by the model using K g 
from Wanninkhof and McGillis [1999] is -38.1 gC/m 2 /yr. In view of this wide range of calculated 
values, we decided to conduct a series of sensitivity runs with the model to assess the effects of differ- 
ent formulations for K g . 

We conducted 7 sensitivity runs: 

(1) baseline run with the best parameter set for the biology components (shown in Table 4) and K o 
from Wanninkhof [1992]; 

(2) same run as in (1) except for using K o from Liss and Merlivat [1986]; 

(3) same run as in (1) except for doubling K o \ 

(4) same run as in (1) except for quadrupling K g ; 

(5) same run as in (4) except for no iron; 

(6) same run as in (4) except for removal of biological drawdown of C0 2 (abiotic); and, 

(7) same run as in (1) except for K g from Wanninkhof and McGillis [1999]. 

The results are summarized in Table 7. A comparison between the sensitivity runs for K (runs 1 , 2, 
3, 4, and 7) reveals that, as K increases, ApC0 2 decreases accordingly but FC0 2 increases much more 
slowly. For example, the mean values for K, ApC0 2 , and FC0 2 from run 1 are 0.34, -20.4, and -38.1. 
Conversely, the mean values for K, ApC0 2 , and FC0 2 from run 7 are 1.287, -2.3, and -45.6, respec- 
tively. Thus, for run 7, K is about 4 times larger, ApC0 2 is 9 times smaller, but FC0 2 is only 20% larger. 
T his means that there is a dynamic negative feedback in the model that prevents the surface flux from 
assuming unrealistically large values no matter how large the value of K is. Table 4 also shows that the 
iron and the biological uptake affect the surface flux significantly. The run without iron limitation (run 
5) provides a surface flux 30% larger, whereas the abiotic run (run 6) reverses the flux direction to 
+15.8 gC/m 2 /yr. The question that still remains is what formulation for K provides the most realistic 
pC0 2 variability. The relative influence of light and nutrient limitation on the phytoplankton growth 
simulated by the model is shown in the climatological vertical profiles of light limitation (L lim ), nitrate 
and ammonium limitation (N ]im ), and iron limitation (Fe Um ) of Figure 11. Within the top 50 meters iron 
is the limiting factor, while below 50 meters, and within the euphotic zone, light is the limiting factor. 

Figure 12 shows the seasonal variations of SST, SSS, nitrate, chlorophyll, total carbon dioxide, in 
situ pC0 2 , and pC0 2 normalized to 10 °C obtained from the model simulation (solid lines) and from 
observations (black circles and triangles). The chlorophyll data consist of two sources. The Wong and 
Chan [1991] data for the period of 1 973- 1978 (triangles), and the National Oceanic Data Center (NOAA / 
NODC) data set for the period of 1959-1980 (circles). All parameters simulated by the model are in 
good agreement with observations, except for the pC0 2 . 
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Table 6. Comparison between model and observed parameters 



Model 

Data 

PP ( eC/m 2 /vr) 

168 

140 

Chlorophyll (mg/m 3 ) 

0.33 

0.4 1 1 
0.32" 

Nitrate (uM) 

10.8 

10.6 

SST (°C) 

7.8 

8.2 

SSS (psu) 

32.68 

32.67 ! 

TC0 2 (mmol/m 3 ) 

2016.1 

2059.7 

Air pC0 2 (fiatm) 

320.9 

329.8 

Ocean pCO, (|iatm) 

300.5 

317.0 

ApC0 2 (patm) 

-20.4 

-12.8 

K (mmol/m 2 /patm/d) 

0.340 

0.1 37 m 

FC0 2 (gC/m 2 /yr) 

-38.1 

-8.3 ni 

-19.1 


I. From Wong and Chan [1991] 

II. From McClain et al. [1996] 

HI. Calculated quantities based on wind speed, C0 2 solubility, and ApC0 2 . 
Values using K from Liss & Merlivat and Wanninkhof, respectively. 



Figure 11. Climatological profiles of light and nutrient limitation predicted by the model. 
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Table 7. Comparison between 7 different model runs. All values were obtained by averaging the 
model-data overlapping years (1973-1978), except for the last run (K from Wanninkhof and McGillis, 
1999) where 21 years (1960-1980) were used for the averages. The first 2 spin-up years (1958-1959) 

were eliminated from the averages 


Run Type 

PP 

(gQntf/yr) 

Chl-a 

(mg/m 3 ) 

no 3 

(M-M) 

Ocean 
pC0 2 Qiatm' 

tco 2 

(mmol/m 3 ) 

ApC0 2 

(patm) 

K=K a 

O 

(mmol/m 2 

/patm/d) 

fco 2 

(gC/m 2 /yr) 

Baseline 1 

168 

0.33 

10.8 

300.5 

2016.1 

-20.4 

0.340 

-38.1 

Baseline2 

168 

0.33 

10.8 

282.3 

2005.4 

-38.6 

0.156 

-29.6 

2xK 0 

168 

0.33 

10.8 

311.4 

2022.5 

-9.5 

0.679 

-42.7 

4xK 0 

00 

\D 

0.33 

10.8 

317.0 

2025.0 

-3.9 

1.358 

-44.4 

No iron 

198 

0.35 

4.6 

314.6 

2023.8 

-6.3 

1.358 

-57.3 

Abiotic 

0 

0 

38.8 

329.7 

2028.3 

7.7 

1.385 

15.8 

K (W 3 ) 

174 

0.33 

10.8 

318.6 

2023.1 

-2.3 

1.287 

-45.6 
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Figure 12. Model (solid line) versus observed (black circles and triangles) seasonal variability 
of sea surface temperature and salinity, nitrate and chlorophyll averaged over the upper 80 
meters, surface total carbon dioxide, pC0 2 at in situ temperature, and pC0 2 normalized to 
10 °C. The gas exchange coefficient of Wanninkhof [1992] was used in this simulation. 
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The model pC0 2 is up to 40 jiatm lower in the winter-spring months and about 20 (latm higher in 
the summer. A significant improvement is achieved when the McGillis and Wanninkhof[ 1999] formu- 
lation for K is used. Figure 13 shows the improved results. Note that the pC0 2 seasonal variability is 
now much closer to the observed data. Figure 14 compares the model and observed pC0 2 interannual 
variability for the period of 1973-1974. Predictions from run 1 (dotted line) and run 7 (solid line) are 
shown together with the observed values (solid black circles). It is quite evident that more realistic 
amplitude and phase of the seasonal cycle are achieved when the McGillis and Wanninkhof [1999] 
formulation for K is used. Therefore we adopted this K formulation for our interannual run. The results 
of the interannual run are discussed in the following sections. 

5.0 SEASONAL VARIABILITY 

Figures 15, 16, and 17 show the seasonal variability (from 1960-1980 monthly averages) of the 
physical and biogeochemical model parameters. These are the vertical velocity, vertical eddy diffusivity, 
temperature, photosynthetically available radiation (PAR), nitrate, ammonium, zooplankton, phytoplank- 
ton, total carbon dioxide, oxygen, and iron. The vertical velocity is very weak in general (maximum of 
3 cm/d) with upwelling peaks in the spring and fall when the wind curl is largest. Maximum downwelling 
occurs in December. The vertical advection effects are minimal when compared to vertical diffusion 
(K ). Maximum surface K v values ranging from 900 to 1000 m 2 /d occur in late fall and winter. The 
depth penetration of large K v values follows the seasonal changes of the mixed layer depth; largest in 
winter (120 m) and smallest in summer (20-40 m). Downwelling light intensity and penetration are 
larger during April-August, peaking in May- June. The wanner temperatures are confined to the top 50 
meters and the time period of May-October, with a peak in August-September. The seasonal variability 
of phytoplankton, nitrate and iron are strongly correlated with top 50-meter peaks extending from May 
through October. This is consistent with photosynthetic consumption of nutrients during the high growth 
season. Total carbon dioxide and oxygen are strongly correlated with temperature due to the strong 
dependence of solubility on temperature. 

Oxygen anomalies relative to the temperature-determined saturation value (Figure 18) show that 
there is a seasonal cycle of air-sea flux, with ingassing in winter and outgassing in summer. The oxygen 
anomaly in the abiotic case is consistently less than in the full coupled model, suggesting that net 
community production ( N p ) is positive at all times of the year. The biologically generated anomalies 
contribute very little to the air-sea flux of oxygen. The difference between the biotic and abiotic cases 
is greatest in summer, with a maximum of 3.3 pM in July, reflecting the seasonal cycle of N p . This 
summer maximum in the NetCP was hypothesized by Wong and Chan [1991] to explain the absence of 
a strong seasonal supersaturation of pCO r The oxygen supersaturation in the summer months is 1.9- 
2.8 times the value expected from thermal forcing alone. This is consistent with the oxygen supersatu- 
rations determined by Emerson et al. [1993], which range from 1.4-2.6 times the values for argon (an 
inert tracer of abiotic effects on oxygen), so rates of biological new production in the model are reason- 
able and even on the high side. 
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Ocean pC0 2 in situ 



1973 1974 1975 1976 1977 1978 1979 


Ocean pC0 2 at 1 0°C 



1973 1974 1975 1976 1977 1978 1979 


Figure 14. Comparison between model and observed pC0 2 (in situ and normalized to 10 °C) 
for the period of 1973-1978. Results from two model runs are shown: (1) the dotted line shows 
the model results using the gas exchange formulation of Wanninkhof [1992]; 

(2) the solid line shows the model results using the gas exchange 
formulation of Wanninkhof and McGiUis [1999], The solid black 
circles are the monthly averaged observed values. 
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Figure 17. Seasonal variability of temperature, total carbon dioxide, oxygen, and iron 

concentrations simulated by the model. 






6.0 INTERANNUAL VARIABILITY 


The interannual variability of the major parameters simulated by the model is illustrated in the 
profile time series shown in Figures 19 and 20, and in the surface time series shown in Figure 21. The 
temperature profile shows distinct warm and cold periods. Two warm periods occur on the series, one 
during 1960-1965 and another during 1976-1980. The period during 1966-1975 exhibits colder (~ 1 °C) 
temperatures in the top 150 meters. The oxygen profile series shows higher (~20 mmol/kg) concentra- 
tions during the cold period due to increased solubility. The phytoplankton concentration in the upper 
20 meters was higher (~0.1 mg/m 3 ) during the warmest temperature events (1961, 1963, and 1979). 

0 2 Anomaly (0 2 -0 2 sat ) 



Figure 18. Seasonal variability of surface oxygen anomaly predicted by the model. 
Results from two runs are presented, an abiotic run (dotted line) and a biotic run (solid line). 


The secular pC0 2 upward trend shown in the top tier of Figure 21 is manifested in the profile series of 
total carbon dioxide concentration in Figure 19. The carbon dioxide increase of about 20 mmol/m 3 in 
20 years, caused by the increased surface flux, is noticeable at depths of up to 250 meters. The highest 
drawdown in the iron, nitrate, and carbon dioxide concentrations occurred during 1961-1963, a period 
of warmest euphotic zone temperatures during which solubility was lowest and phytoplankton growth 
was highest. Figure 22 shows the yearly and depth averaged (0-100 m) temperature and TC0 2 variabil- 
ity, and the surface C0 2 flux and oceanic pC0 2 variability for the period of 1960-1980. The interannual 
changes and trends are highlighted in this figure. Note the cooling trend in the temperature during the 
1960s, followed by a warming trend during the mid to late 1970s. The air-sea C0 2 flux is always 
negative, indicating that OWS P is a sink of atmospheric CO r The pC0 2 trends upward almost mono- 
tonically showing an increase of about 30 (latm in 20 years. The TC0 2 follows the pC0 2 trend, show- 
ing that the atmosphere and ocean are tightly coupled at OWS P. 


26 




Depth (m) Depth (m) Depth (m) Depth (m) 






Depth (m) Depth (m) Depth (m) Depth (m) 




60 62 64 66 68 70 72 


76 78 80 


Figure 20. Interannual variability of phytoplankton, iron, nitrate, and ammonium 

simulated by the model. 
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Figure 21. Time series of simulated atmospheric pC0 2 , ocean pC0 2 , SST, mixed layer depth, 
total carbon dioxide, oxygen, iron, nitrate, and ammonium. 
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7.0 CARBON FLUX BUDGET 

Figure 23 represents the climatological (1960-1980) ecosystem carbon flux balance. The main 
compartments in this flow chart are the phytoplankton, zooplankton, respired DOC, and total C0 2 
carbon stocks. The numbers indicate the carbon flux between ecosystem components in gC/m 2 /yr. The 
ammonium and nitrate boxes were replaced by the net community production (76 gC/m 2 /yr), which 
represents all sources and sinks of macronutrients combined (equation 14). Note that the air-sea (44 
gC/m 2 /yr) and bottom (23 gC/m 2 /yr) carbon fluxes are balanced by the POC+DOC export (60 gC/m 2 / 
yr) and fecal pellet loss (7 gC/m 2 /yr) through the bottom. The total gross production (uptake) is 174 gC/ 
m 2 /yr. Therefore, the total carbon required for photosynthesis is partitioned among the following sources: 
44% (76 gC/m 2 /yr) originates from the TC0 2 pool - of these 44%, 5% comes from DOC respiration, 
13% comes from the bottom flux, and 25% from the atmosphere (air-sea flux). The other 56% (98 gC/ 
m 2 /yr) of the required equivalent carbon comes from the recycling of zooplankton and phytoplankton 
losses (death, respiration, and fecal pellet production). 

8.0 SUMMARY AND CONCLUSIONS 

A coupled ecosystem/carbon flux model was developed, validated, and used to simulate bio- 
geochemical parameters at OWS P. A series of sensitivity runs revealed that the most significant 
improvement in the simulation of the surface pC0 2 is obtained by the use of a recently developed cubic 
formulation for the gas exchange coefficient (Wanninkhof and McGillis, 1999). This result shows that 
the ocean and atmosphere are much more closely coupled at high (> 10 m/s) wind speeds than indi- 
cated by the previous wind-dependent formulations (Liss and Merlivat, 1986; Tans et al., 1990; 
Wanninkhof, 1992). 

All biogeochemical parameters, when averaged over the period of available concurrent observa- 
tions (1973-1978), are within less than 5% of the observed values. The only exception is the sea-air 
C0 2 flux which departs more significantly from previously reported values. We attribute this differ- 
ence to the elusive nature of sea-air C0 2 estimates in the ocean and the variety of formulations used in 
the literature to derive the gas exchange coefficient. For example, our lowest C0 2 flux estimate is -29.6 
gC/m 2 /yr when we use the Liss and Merlivat [1986] gas exchange formulation. Wong and Chan [1991], 
also using the Liss and Merlivat [1986] gas exchange formulation, report a value of -8.4 gC/m 2 /yr. Our 
C0 2 flux value is even higher (-45.6 gC/m 2 /yr) when we use the Wanninkhof and McGillis [1999] gas 
exchange formulation. However, the model estimates for the surface flux and those derived from ApC0 2 
measurements are not entirely equivalent because the model provides a dynamically coupled flux esti- 
mate, whereas the air-sea flux estimates from observed ApC0 2 are obtained from a relatively small 
number of discrete data points. Therefore, global estimates of air-sea C0 2 flux derived using air-sea 
flux models to extrapolate from observed pC0 2 may not be realistic. In the short term, fully prognostic 
models may give the best estimates provided that they achieve sufficient agreement with other observ- 
able quantities such as temperature, chlorophyll, and TC0 2 . 
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Figure 22. Yearly averaged time series of ocean temperature (0-100 m mean), air-sea carbon 
dioxide flux, surface carbon dioxide partial pressure, and total carbon dioxide concentration 
temperature (0-100 m mean) simulated by the model for the period of 1960-1980. 










CLIMATOLOGICAL ECOSYSTEM CARBON FLUX BALANCE 
(g C/m 2 /yr) 

Figure 23. Flow chart of simulated climatological carbon flux balance showing the principal 
carbon exchange compartments in the coupled ecosystem/carbon-flux model. 
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